#ifndef __GLEW_H__
#include <GL/glew.h>
#include <GL/wglew.h>
#endif
#include "StimDotGroundDeperspective.h"

#include <iostream>

CStimDotGroundDeperspective::CStimDotGroundDeperspective(const bool& bAlsoCeiling_, const enum Stims::ObjType& eObjType_/* = Stims::GLpoint*/)
    :
CStimDotGround(bAlsoCeiling_, eObjType_) // set dot volume type and dot type
{
}

CStimDotGroundDeperspective::~CStimDotGroundDeperspective(void)
{
}


std::string CStimDotGroundDeperspective::GetName()
{
    if (_eStimType==Stims::GroundCeiling)
        return std::string("Ground & Ceiling, deperspective");
    else
        return std::string("Ground, deperspective");
}

void CStimDotGroundDeperspective::InitVectors()
{
    // always
    _vvPointList    .clear();
    _vvPointList    .resize(_nNArea+1);
    _vNumDot        .clear();
    _vNumDot        .resize(_nNArea+1);
    _vHalfLengthX   .clear();
    _vHalfLengthX   .resize(_nNArea+1);
    _vHalfLengthY   .clear();
    _vHalfLengthY   .resize(_nNArea+1);
    _vHalfLengthZ   .clear();
    _vHalfLengthZ   .resize(_nNArea+1);

    // limited lifetime
    _vvPointLifeList.clear();
    _vvPointLifeList.resize(_nNArea+1);

    // when saving dots to file, vector of point IDs
#if STIM_DOT_SAVE_DOTS
    _vvPointIDs.clear();
    _vvPointIDs.resize(_vvPointList.size());
#endif
}

void CStimDotGroundDeperspective::ComputeDistribution()
{
    // initialize
    std::vector<double> vDensity(_nNArea+1,0.);

    // compute frustum definition
    if (_dNearPlane < _dEyeHeight/_dTanFOVy_2)
    {
        // correct near clipping plane to what is visible inside screen
        _dNearPlane = _dEyeHeight/_dTanFOVy_2;
        std::cout << "Near Plane adjusted to: " << _dNearPlane << "m" << std::endl;
    }
    double dHighVertical    = _dEyeHeight/_dFarPlane*_dNearPlane;           // part of imaging plane behind far clipping plane

    // To create our stimulus, we need to keep the density on the imaging plane equal.
    // So, first compute the vertical extend of the imaging plane. As it is put at
    // the near depth plane, it always ranges from 0 to HighVertical, the maximum height
    // on imaging plane. We then divide this range into equal steps.
    double dStep  = (_dEyeHeight-dHighVertical)/_nNArea;

    // each rectangle represents an equal area on the imaging plane, thus divide total number
    // of visible dots by the number of rectangles used
    double dNumPointPerImagingRegion = _nNPoint/double(_nNArea);

    for(int i=0; i<_nNArea+1; i++)
    {        
        _vHalfLengthZ[i] = (i*dStep)*_dNearPlane/(_dEyeHeight - i*dStep) + _dNearPlane;        // backproject point on imaging plane to get length from origin to end of rectangle (as rectangle are centered on the observer, thats half the length of each dimension of the rectangle)
        _vHalfLengthY[i] = _vHalfLengthZ[i]*_dTanFOVy_2;                            // knowing FOV and aspect ratio, we can compute the Y also, Y needed as needed when testing if point inside frustum
        _vHalfLengthX[i] = _vHalfLengthZ[i]*_dTanFOVx_2;                            // knowing FOV and aspect ratio, we can compute the X also
        if(i==0)                                                                    // density inside the smallest rectangle is zero, as it covers the area before the near clipping plane
            continue;

        double dSmallArea = _vHalfLengthZ[i-1]*_vHalfLengthX[i-1];                  // Area of previous rectangle
        double dLargeArea = _vHalfLengthZ[i]*_vHalfLengthX[i];                      // Area of this     rectangle
        double dAreaTrapezoid  = dLargeArea - dSmallArea;                           // Area exclusive to this rectangle

        vDensity[i] = dNumPointPerImagingRegion/dAreaTrapezoid;                     // corresponding density
    }


    // key point of the algorithm (rephrase the below once we've figured it out!):
    // Density now also include next rectangle, so subtract density to make density correspond to number of dots that should be in this rectangle
    for(unsigned int i=1; i<vDensity.size()-1; i++)
    {
        vDensity[i] -= vDensity[i+1]; 
    }

    // given density and area, we can compute how many dots each square should have
    double _dNearTriangleArea = _dNearPlane * _dNearPlane*_dTanFOVx_2;                        // area of the triangle behind the near clipping plane (near observer) -- halfz*halfx == halfz*x / 2
    for(int i=0; i<_nNArea+1; i++)
        _vNumDot[i] = int((_vHalfLengthX[i]*_vHalfLengthZ[i] - _dNearTriangleArea) * vDensity[i] + .5);   // how many dots inside this trapezium

    // number of dots in stimulus
    int nCount = 0;
    for(int i=0; i<_nNArea; i++)
        nCount += _vNumDot[i];

    std::cout << nCount << " dots in stimulus" << std::endl;
}